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Abstract 

We present a general scheme for the computation of the time dependent 
(TD) quadratic susceptibility (x^) of an extended insulator obtained by ap- 
^ ' plying the '2n + 1' theorem to the action functional as defined in TD density 

functional theory. The resulting expression for x^ includes self-consistent 
local-field effects, and is a simple function of the linear response of the sys- 
tem. We compute the static % °f nine III-V and five II- VI semiconductors 
using the local density approximation (LDA) obtaining good agreement with 
experiment. For GaP we also evaluate the TD x^ for second harmonic gen- 
eration using TD-LDA. 
42.65.Ky,71.10.+x,78.20.Wc 
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Nonlinear optics is a growing field of research which has applications in many technical ar- 
eas such as optoelectronics, laser science, optical signal processing and optical computing [|TJ. 
In these fields the description of several physical phenomena, such as optical rectification, 
wave-mixing, Kerr effect or multi-photons absorbtion, relies on the knowledge of the non- 
linear optical (NLO) susceptibilities. Moreover nonlinear spectroscopy is a powerful tool to 
analyze the structural and electronic properties of extended and low dimensional systems. 
In the present work we give a general scheme to compute from first principles the time de- 
pendent (TD) quadratic susceptibility (x^) of real materials within TD-density functional 
theory (DFT). Futhermore we show that the values of the static x^ obtained in the local 
density approximation (LDA) are in good agreement with measured values for the cubic 
semiconductors. Our approach makes feasible the computation of x m ce hs containing up 
to an hundred atoms, since it requires the same numerical effort as the computation of the 
total energy. This allows the evaluation of x^ f° r systems of technological and scientific 
relevance which can not be handled by the traditional methods, such as surfaces or crystals 
of organic molecules. 

Nowadays many first-principle calculations for the ground state properties of materials 
are performed within DFT. Even in its simplest form, namely in the LDA for the exchange 
and correlation energy this scheme gives results which, in many cases, are in surprisingly 
good agreement with experiments. A rigorous extension of DFT to TD phenomena has 
been proposed in Ref.s @,[|]. Although the available approximations for the exchange and 
correlation energy are less accurate in the TD domain than in the static case, this scheme is 
sufficiently general to allow many possible improvements in the future. Therefore TD-DFT 
seems to be a promising framework for the study of the NLO susceptibilities. 

Standard quantum-mechanical perturbation theory can be used to compute the x^ 1 ■ 
The straightforward application of perturbation theory leads to an expression for x^ , which 
diverges for an infinite solid in the static limit. However, for an insulator, these divergences 
have been shown to be apparent . This kind of approach has been applied to compute the 
■^(2) from first principles. The non self-consistent expression for x reported in Ref. @] has 
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been evaluated by Huang and Ching || using the DFT-LDA wavefunctions and eigenvalues. 
A fully self-consistent theory of the NLO susceptibility within DFT has been proposed 
in a series of papers by Levine and Allan ||. Their method is feasible but algebraically 
very involved due to the necessity of dealing with the second order perturbation of the 
wavefunctions and with the apparent divergences. Their final expression is not easy to 
handle and its evaluation requires summations over the conduction band states, which are 
time consuming and difficult to converge. 

In a previous paper two of us |7j have shown that it is convenient to regard the static 
as a third order derivative of the total energy with respect to an uniform electric field. 
We pointed out that this derivative can be obtained by combining a Wannier representation 
of the electronic wavefunctions with the '2n + 1' theorem of perturbation theory [§|§. We 
also found an equivalent expression of the static x m terms of Bloch wavefunctions. 

In the present letter we show that the method of Ref. [[/[ applies also to TD periodic 
perturbations and to the self-consistent TD-DFT functional. The TD x^ can be regarded 
as a third order derivative of the total action. The stationary principle for the action func- 
tional |2|||], which replaces in the TD case the miminum principle for the energy functional, 
allows the use of the '2n + 1' theorem. As in the static case the third order derivative de- 
pends only on the unperturbed wavefunctions and on their first order change due to the TD 
electric field. All the self-consistent contributions are included in the formalism in a simple 
way. The final expression avoids perturbation sums and does not present any apparent di- 
vergency. We apply our formalism to the computation of the the static x^ °f nme HI-V 
and five II- VI cubic semiconductors within the LDA. For GaP we also evaluate the TD x 
for second-harmonic generation (SHG) using TD-LDA 1^ . 



In the Kohn and Sham (KS) formulation of DFT the ground state density n 9 (r) of a 
system of N interacting electrons in an external potential V ext (r) is written in terms of N/2 
single particle wave-functions {4> 9 }- The set {4> 9 } minimizes the KS energy functional 
and the ground state energy is obtained as E 9 = E[{(p 9 }]. A formalism similar to that of the 
static case can be introduced also in the TD domain if one restricts to Hamiltonians periodic 



in time and to the evolution of the system which is steady and has the same periodicity of the 



Hamiltonian [|TTJ. In TD-DFT the TD steady density n s (r, t) of a system of N interacting 
electrons in an external TD potential V ex t(r, t), periodic in time with period T, is expressed 
in terms of a N/2 TD single particle wave-functions {ijj s } @||- The set {ip s } make stationary 
the KS action functional A[{^}], i.e. 

6A[W}]/6(Mt)\ = 0> (!) 
and the steady action is obtained as A s = A[{tp s }]. The KS action functional is 

defined as (atomic units are used throughout) : 

r T rf+ \ N I 2 i a r 

A[{1>}) = 1 J E2(^(t)|--V 2 -^|^(t))+yrf 3 r^(r,tMr,t) 

+A H [n] +A xc [n\. 

Here \^(t)) = Mt + T)), (^(t)l^-(t)) = Sy, A H [n] = / T dt/T J d 3 rd 3 r'n(r, t)n(r', t)/(2\r - 
r'|) is the Hartree functional, A^n] is the exchange and correlation functional, and 
n(r, t) = ThIx 2(^(t) |r) (r |^(^)) , where the 2 factor is for spin degeneracy. At this stage no 
approximation for the exchange and correlation functional is made. The stationary principle 
in Eq. (|l|) yields the TD KS equations: 

here e k are the steady states eigenvalues, H KS (t) = — |V 2 + V ext (r,t) + VHxc[n}(r, t) is the 
time dependent KS Hamiltonian, and V^d/i] (r, t) = T5(A H [n] + A xc [n])/Sn(r,t). 

Now we consider a potential of the form V ext (r, t, a) = V e ° t (r) + a\ei ■ r cos(co>it) + 0262 • 
rcos(a;2t) + 0363 ■ rcos^t), where ei, e 2 , e 3 are unit vectors describing the orientation of 
three TD uniform electric fields, u\ + u> 2 + uj 3 = and a = (ai, a 2 , a 3 ) describes the strength 
of the fields. Then the steady state wavefunctions {^"(a)} and action A s (&) depend also on 
a. Note that for a = the potential is time independent and the action coincides with the 
static DFT energy. By using the Hellmann-Feynman theorem we obtain the derivative of 
the action with respect to the parameter a\. 
dA s (z) f T dt 



f dt f 

/ — cos(a;it) / d 3 r e 1 ■ r n s (r, t, a) = — ei • P s (co>i, a) V, 
Jo T J 



where V is the volume of the system and P s (u;i, a) is the macroscopic electronic polarization 
per unit volume, oscillating at frequency U\ fl^. Then the quadratic susceptibility tensor, 



which is defined as x2 ) ; e 2 ,e 3 ( _a; i; ^2, w 3 ) = y^S^f 1 , is equal to 



(2) 



2 d 3 A s (0) 



;e2,e3 



V da\da2da^ 

The computation of the derivatives of A s (&) with respect to a, can be performed by 
using the '2n + 1' theorem which states that the derivatives up to order 2n+ 1 of the steady 
action depends only on the change of the orbitals up to order n: 



(2) 



Q a 2n+1 y Q a Q a n J ' 

where V 2n+1 is a polynomial of degree 2n + 1 in its arguments. Indeed, as shown in [|7].^| . 
Eq. (H) relies just on the stationary condition, Eq. (p. Therefore Xeve 2 e 3 ( — u 2i = 
— yV 3 ( 8 ^a & ^ ) • The derivation of an explicit expression of V 3 for an infinite periodic 
system requires a particular care because the expectation value of the r operator between 
Bloch states is ill defined. In an insulating solid this problem can be solved following Ref. [|7|]: 
first we apply the '2n + 1' theorem in a Wannier representation where the r operator is 
well defined, then we recast the resulting expression in a Bloch representation. The final 
expression is: 



N/2 



m,n (7=± 



d 3 k 



BZ (27r) ; 



( M k,m 



e 2 



-id 



0k y U k,n) i U lc,n \) \ U k,m) 

+fim,n(Uk,n \^Hxc\ U k,rn) ~ ( U k,m I ^Hxc\ U k,n) ( U k]n \ U k,m) 
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J d 3 rd 3 r'd 3 r"K xc (uj 2 , u 3 , r, r', r")n ai (r)n a2 (r> a3 (r") 



+n{i,2,3}. 

Here II{1, 2, 3} indicates the sum over the 5 permutations of the indexes 1, 2, 3, and 
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(3) 



V%Jt) = J d 3 r> 

N/2 

™ ai (r) = 2£E 

m <r=± 



r — r 

Qd 3 k 
bz (2tt) 3 



- + M cc (w 1 ,r,r / ) 



n ai (r') 



Re 



(<ml r )( r Km) 



Jo dnlr, onir, t) 



K xc {u 2 ,Wz,r,Y',r") = f 

Jo 



5n(r,Q)5n(r',t) 
t ^ 5M xc [n°}(u 2 ,r,r>) 
8n(r",t) 



n° is the unperturbed charge density, \u^ m ) is the periodic part of the unperturbed Bloch 
eigenstate normalized on the unit cell Q, with eigenvalues e£ m , and \u^!^) are the perturbed 
orbitals projected on the unperturbed conduction band subspace, i.e. the solution of the 
linear system: 

(elm - His ± ^l)K^ ) = + V Hx)j K, m )- (4) 



with Q k = l-E! 2 |< m )K 



k,m I 



Note that the evaluation of Eq. (|3|) requires only the knowledge of unperturbed valence 
wavefunctions \u^ m ) and of their linear variation li^ 1 ^)- Moreover the solution of Eq. 
can be obtained by minimizing a suitably defined functional with a numerical effort similar 
to the computation of the total energy PjI5|,?|. Thus our formulation makes the evaluation 
of m systems containing up to an hundred atoms feasible. 

We have applied Eq. (|3|) to compute the static of nine III-V and five II- VI cubic 
semiconductors, evaluating the exchange and correlation energy within the LDA. We do 
not use any scissor operator to correct for the LDA band-gap error, contrary to what has 
been done in other ab-initio calculations |5]J|. Indeed the static x^ is a ground state 
property, which is defined as a difference of ground state total energies and it is not related 



to the LDA band gap [15]. We think that improvements over LDA require a better E x 



functional, which could be ultra-nonlocal [16], instead of an ad-hoc correction of the LDA 



band-gap. Furthermore our purpose here is to give reference values for the static x^ which 
are completely consistent within LDA. 

We used norm-conserving pseudopotentials and a plane-wave kinetic energy cut-off of 
24 Ry. The derivative with respect to k which appears in Eq. (|3]) has been computed by 
means of finite differences. We have found that the effect of d electrons is important for Ga 
and In atoms, and it necessary at least to use the nonlinear core corrections (NLCC) |20| to 
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obtain the correct LDA values for x^ m the compounds containing these elements pTJ . For 



II- VI semiconductors the effect of the cation d electrons is even more important p2| and our 
reported values have been computed using the NLCC. For A1P, AlAs, GaP and GaAs we 
have also verified that our results for the x^ reproduce the LDA values obtained in Ref. 
if the same pseudopotentials (without NLCC) and lattice constants are used. 

In Table | we report the values of the x °f the III-V and II- VI cubic semiconductors 
computed at the theoretical LDA lattice constant (a ), also reported in the Table. On the 
same Table we show also the direct band-gap at the T point, E^, and the static dielectric 
constant e^. Known experimental values for ao, and 6^ are reported in parenthesis. 
Well established experimental data for x^ do not exist since the values reported by different 
authors may differ by more than a factor of 2. Moreover, in some cases only data obtained 
at frequencies close to the absorbtion edge are available. Therefore we refer the readers to 
Ref.s [0,|19|,§] for a complete review of the experimental results. Just to give an indicative 
value, we show in parenthesis the experimental results from Ref. [plj which correspond to the 
lower frequencies. For GaP, GaAs and CdSe we have taken the values from Ref. |19| obtained 
after an appropriate rescaling of the experimental data. In the case of InAs we cannot 
compute x an d £oc since within LDA the system is a metal. For all other compounds the 
computed x^ are i n the range of variation of the available experimental data [[0|,[n^,[J . 

As a second application we compute the TD x^ for SHG of GaP. For this calculation we 
used the TD-LDA. In the TD case the use of LDA is less justified since, in general, it does not 
describe correctly the position of discrete excited levels and absorption edges as difference 
to the exact TD-DFT We note that this is a limitation of the approximation to A xc 
used here, and not of Eq. (§) itself. Since LDA is not expected to perform sufficiently well 
in the TD domain we have used the pseudopotential without NLCC which at its theoretical 
lattice constant (a = 10.01 a.u.) gives a gap (Er = 2.8 eV) and thus an absorption edge, 
which is incidentally close to the experimental one. 

In Table |I| we report x^ 2 ^(2c(j; cu, u;) computed as a function of lo in the non-absorbing 
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regime. The experimental measurements are taken from Ref.s [p!9|,^1 . 

In conclusion we have presented a consistent theory for the computation of the static 
and dynamic nonlinear optical susceptibilities within DFT. To this purpose we have ap- 
plied for the first time the '2n + 1' theorem to the TD-DFT action functional. We have 
presented applications to cubic semiconductors. Our results show that LDA reproduces the 
experimental static nonlinear susceptibilities in these compounds without using any scissor 
operator, provided that the computations are performed at the theoretical lattice constant 
and NLCC are included for Ga, In, Zn and Cd atoms. 
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TABLES 

TABLE I. LDA nonlinear susceptibilities (x^) of IILV and II- VI cubic semiconductors. We 
report also the theoretical lattice constant (ao), the direct gap at the L point (Er), and the dielectric 
constant (£<x>)- Experimental values are given in brackets. All computations are performed with 
28 special k-points, but for InSb for which we used 60 special k-points. 
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